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Abstract. We develop a rigorous and highly accurate technique for calculation of 
the Berry phase in systems with a quadratic Hamiltonian within the context of the 
Kitaev honeycomb lattice model. The method is based on the recently found solution 
of the model which uses the Jordan- Wigner-type fermionization in an exact effective 
spin-hardcore boson representation. We specifically simulate the braiding of two 
non-Abelian vortices (anyons) in a four vortex system characterized by a two-fold 
degenerate ground state. The result of the braiding is the non-Abelian Berry matrix 
which is in excellent agreement with the predictions of the effective field theory. The 
most precise results of our simulation are characterized by an error on the order of 10~^ 
or lower. We observe exponential decay of the error with the distance between vortices, 
studied in the range from one to nine plaquettes. We also study its correlation with 
the involved energy gaps and provide preliminary analysis of the relevant adiabaticity 
conditions. The work allows to investigate the Berry phase in other lattice models 
including the Yao-Kivelson model and particularly the square-octagon model. It also 
opens the possibility of studying the Berry phase under non-adiabatic and other effects 
which may constitute important sources of errors in topological quantum computation. 
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1. Introduction 

Quantum statistics is intimately related to the adiabatic exchange of indistinguishable 
particles. Exchanging two particles twice results in a loop trajectory which in three 
dimensional space can be smoothly contracted to a point, equivalent to no trajectory. 
The particles' wavefunction thus remains unchanged after two subsequent exchanges, 
and after one exchange can transform either in a symmetric or an antisymmetric fashion, 
giving Bose-Einstein and Fermi-Dirac statistics respectively. In two dimensional space 
such a contraction of the particles' trajectory is impossible. This gives rise to a different 
type of quantum statistics. Algebraically, adiabatic exchange operators correspond to 
the elements of the permutation group Sn for three and higher dimensional systems, 
and to the elements of the braid group -B^v for two dimensional systems [1,2]. Both 
groups are formed by — 1 generators ri, r^^i, obeying the constraints 



where the generator interchanges the two particles at the positions i and i + 1. The 
quantum statistics then arises from the unitary irreducible representations (irreps) of 
these groups [1]. The group possesses two one-dimensional irreps which correspond 
to bosonic and fermionic statistics. Its higher dimensional irreps can be replaced by 
bosonic and fermionic statistics when a hidden degree of freedom is introduced [3]. 

On the other hand, one- dimensional representations of the braid group Bn can be 
any phase e*^ where 6 G [0,27r) (hence the name anyonic) [4,5]. More interestingly, the 
braid group also permits multi-dimensional unitary irreducible representations which 
give rise to non-Abelian statistics. Any exchange of particles then leads to a unitary 
rotation of a state vector of the system within a D-fold degenerate ground state. The 
degree of degeneracy D depends only on the presence of A^ well-separated identical 
particles. As this is solely linked to the topology of the underlying configuration space, 
braiding the particles is the only way to induce nontrivial operations within this ground 
state subspace. Consequently, the system is immune to any local perturbations or 
fluctuations as long as these do not exceed the spectral gap which separates the rest 
of the system's spectrum from this decoherence free subspace [6]. This capability to 
implement unitary operations within an intrinsically fault-tolerant framework offers 
promising applications in quantum information processing, specifically topological 
quantum computation [6,7]. 

A larger body of research results suggests that non-Abelian anyons are physically 
realized as localized quasiparticle excitations of many-body systems. These for example 
include the systems that manifest the fractional quantum Hall (FQH) effect [8, 9], 
Px + iPy superconductors [10] and the Kitaev honeycomb spin lattice model [11] (with 
its proposed realizations [12,13]). Theoretical studies show that all these systems can 
be effectively described by a similar topological quantum field theory (the Ising and 
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the related SU {2)2 Chern-Simons theory) which can be characterized by three particle 
types labeled as 1 (i.e. vacuum or trivial topological charge), e (fermion) and a (the 
non-Abelian anyon). These satisfy certain fusion and braiding rules which will be 
specified later. The experimental study of non-Abelian anyons is indeed of fundamental 
importance and so far the experimental observations have yielded encouraging results, 
though it is to be said that final verification of non-Abelian statistics remains a great 
challenge [14]. 

Physically, anyonic statistics arise from the evolution of the system under adiabatic 
interchanges of these quasiparticles. According to the adiabatic approximation, a 
physical system remains in its instantaneous eigenspace if a given perturbation is acting 
on it slowly enough and if there is a gap between the corresponding eigenvalue and 
the rest of the Hamiltonian's spectrum. When these perturbations draw a smooth and 
closed trajectory C(A) in the parameter space, the unitary evolution of the system in 
the n dimensional eigenspace is given by the Berry phase (matrix) B{C) 



where |$'(A)) and |$'^(A)) are the eigenstates of the system's Hamiltonian at the value 
of the parameter A. 

In this paper, we directly evaluate the non-Abelian statistics of the Ising anyons of 
the Kitaev honeycomb model. In particular, we numerically calculate the non-Abelian 
Berry phase (matrix) which governs the evolution of the system under the adiabatic 
exchange of two a-particles (vortices) of the Kitaev honeycomb model. This work can 
be seen as an accurate test of the non-Abelian statistics in the Kitaev model which 
offers applications in the context of topological quantum information processing and 
computation. Moreover, it provides a direct way to study the non-Abelian statistics 
in lattice models with a quadratic fermionic Hamiltonian and as such it complements 
similar efforts carried out in the context of continuous systems [15-17]. 

The non-Abelian Berry phase calculations in the Kitaev model have been a subject 
of study by Lahtinen and Pachos [18]. They developed an interesting technique for 
inducing the vortex motion in the Kitaev honeycomb lattice model which we have 
utilized in the present work, though within a different solution of the model. While 
the previous studies established the non-Abelian nature of the statistics, the results on 
both the exact form of the braid matrix and on the exponential convergence with vortex 
separation were not conclusive. 

We establish these results rigorously for much larger vortex separations and to a very 
high accuracy, but we would like to emphasize that our approach goes beyond a mere 
technical improvement. Our calculations rely on the solution of the Kitaev model which 
was presented recently by Kells et al. [19]. This solution employs the Jordan- Wigner- 
type fermionization in the exact effective spin-hardcore boson representation of the 




(4) 
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model and uses no redundant degrees of freedom, thus allowing us to work directly with 
physical eigenstates of the system. This allows us to calculate the Berry phase associated 
with braiding vortices at the minimal distances for up to nine plaquettes. The simulation 
also reaches a very high degree of accuracy as measured by the Frobenius distance 
between the Berry matrix obtained from our simulation and the exact Berry phase from 
the effective field theory. The accuracy of our calculations increases exponentially fast 
with the vortex distance, achieving results which are characterized by errors on the order 
of 10 ~^ and lower. 

We thus present in this work a very accurate technique for the non-Abelian Berry 
phase calculation. The accuracy of the calculations allows us to see the exact dependence 
of the simulated Berry phase on the details of the model, like the splitting of the ground 
state level which is intrinsic to any finite system. It also shows the dependence on an 
exact implementation of the braiding operations, potentially allowing for an analysis 
of non-adiabatic effects, similar to that provided quite recently in a different context 
in [20]. The importance of these effects derives from the fact that they will likely 
be a crucial source of errors for any topological quantum information processing and 
computation. Possible applications thus extend to modeling and implementation of 
quantum information protocols whose reliability can be tested under various effects; 
for example, disorder. Naturally, the first step in this potentially fruitful story is the 
demonstration of highly accurate and sufficiently large scale direct simulation of the 
Berry phase in the Kitaev model, as presented in this manuscript. 

The results we present show strong agreement with the statistical properties of Ising 
anyons derived from the effective theory. One can naturally wonder what the meaning 
of calculating the Berry phase is if we already know it from the relevant effective field 
theory. The point here is that the effective theory gives us nearly no ground to test the 
stability of the Berry phase under the effects mentioned above or to make any conclusions 
about its implementation under (more) realistic conditions and about what pitfalls 
we should expect in such situations. Moreover, the accurate calculations like those 
presented here provide important predictive power in the analysis of topological phases 
in other lattice models including, for example, the Yao-Kivelson model [21] and the 
square-octagon model [22] which exhibits a kaleidoscope of topological phases including 
the Ising and SU{2)2 phases. We believe that this is highly relevant to implementation 
of quantum information processing in Majorana fermion systems [23]. 

The paper is organized as follows. In Section 2, we review the exact solution of 
the Kitaev model on a honeycomb lattice by explicitly describing its eigenstates based 
on the method presented in [19]. We first define Jordan- Wigner types fermions in the 
honeycomb lattice and then represent the original Hamiltonian using these fermions. 
The resulting Hamiltonian is in a quadratic fermionic form which can be solved exactly. 
Then in Section 3, we discuss how to smoothly move the vortices within our solution 
of the model and investigate the adiabaticity of the anyonic motion. When the anyons 
follow a cyclic path adiabatically, the evolution of the system is governed by the Berry 
matrix; the numerical method for calculating this matrix will be presented in Section 
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4. The presented method can be apphed to the Berry matrix of any system having a 
quadratic fermionic Hamiltonian. In Section 5, we discuss the results of the numerical 
calculation by comparing them with the expected statistics and present an error analysis 
of the numerical calculations. 

2. The Kitaev Honeycomb Model 




Figure 1. (Color online) Kitaev honeycomb lattice model (more details in the text). 
Vertex coloring emphasizes the bipartite lattice structure. The operators Wp, Sq^,, 
and Ly are defined as products of single- and two-body terms. 



The Kitaev honeycomb model is a spin-1/2 lattice model in which spins are located 
on the vertices of a honeycomb lattice (see Figure 1(a)), and it has the following 
Hamiltonian: 

x-links {/-links z-links 

where i and j are the position indices of the spins, , a = x,y,z, are the coupling 
coefficients of the two-body interaction operator K°'j = ufa'j on the link (^, j), and the 
erf are the Pauli operators. 

The model is exactly solvable and contains three equivalent gapped A phases for 
parameters satisfying .P > .P + , .P > + or > P + J^, and a gapless B 
phase for the other values of the parameters. Furthermore, adding to the Hamiltonian 
(5) a term which breaks the time-reversal and parity invariance of the model opens a 
spectral gap in the B phase and allows the realization of non-Abelian anyons of the Ising 
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type. This additional term, defined as 

6 

v = j2pl, = Pp+pS+ Pi + Pi + Pi +pi ( 6 ) 

1=1 

,^1 I ,2 I ,^3 I ,A ^x 

+ 4 ^3^I^I + 4 

where p indexes honeycomb plaquettes, represents the time-reversal and parity 
invariance breaking effect of a weak magnetic field [11] 

y=- 5Z(^xCt| + hyo] + (7) 
i 

as it emerges from perturbation theory on the third order. The coefficients of the 
effective term / at the plaquette "p (Figure 1(a)) are related to the magnetic field as 
K ~ h^K^z J ^ jx _ jy _ jz what follows we will consider only the effective 
magnetic field (6). 

The model has a commuting set of plaquette operators 

LT/ — — rr^ 

VVp .— -fV;^ 2-"-2,3-"-3,4-"-4,5-"-5,6-"-6,l — " i O 2 (Tg (Tg (Tg 

for each hexagon p which also commute with the total Hamiltonian Htot = "H + V. The 
eigenvalues of Wp correspond to whether the plaquette p is occupied by a vortex (— 1) 
or not (+1). The vortices carry unpaired Majorana modes for odd values of the Chern 
number u [11]. For translationally invariant configurations of the Kitaev honeycomb 
model, it is found that z/ = ±1 depending on the direction of the magnetic field [11]. 
Majorana modes exhibit non-Abelian statistics [24] corresponding to cr-particles of the 
Ising anyons. We will discuss their properties in more detail in Section 5. 

The model can be solved by various fermionization techniques, but here we will 
use the solution introduced in the paper [19] for our purposes. This solution has the 
advantage of giving the eigenstates of the system explicitly whereas it is not practically 
possible to do that in Kitaev's original solution. The original solution maps the 
spin degrees of freedom of the model to Majorana fermions. This requires each spin 
degree of freedom to be embedded in an extended Hilbert space of four dimensions. 
Obtaining physical states then requires projections from the eigenstates of the extended 
Hamiltonian which is hard to achieve in practice and thus limits the extent of numerical 
calculations. For the sake of self-completeness and clarity for further arguments, we 
start with a brief discussion of the solution we will rely on. 

2.1. The exact solution of the model 

Here we focus on x A^^ lattices on a torus where and Ny are the numbers of 2;-links 
in the and riy directions respectively, as in Figure 1(a), where the dotted lines define a 
4-by-4 lattice whose opposite sides are identified. Let us label the ^-hnks by g = (g^, %) 
with respect to the z-link at the origin O = (1, 1). Spins are denoted by either empty o 
or full • circles, reflecting the bipartite structure of the underlying honeycomb lattice. 
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Periodicity on the torus imposes Yip = 1 on the plaquette operators and gives rise 
to homologically nontrivial loop symmetry operators and Ly (see Figure 1(b), (c)) 
which have ±1 eigenvalues and commute with all Wp and the total Hamiltonian [25]. 
Therefore, for a 2A^-spin system on a torus, there are + 1 independent operators 
which split the total Hilbert space of the system into 2^+^ different 2^~ ^-dimensional 
subspaces. 

Before we define the fermions on the lattice, let us define the string operators 
between an arbitrary location on the lattice q = {qxiQy) and the origin O = (1, 1) (see 
Figure 1(a)) as Sg^, := SySxCr^ , and Sg^o '■= <7g,o^g..'S'g,» where the string denotes the 
successive applications of c^c^ and a^a^ to the 2;-links and x-links of the interval [O, 
{qx, 1)) respectively and similarly Sy is the successive applications of a^a^ and cr^cr^ to 
the 2;-links and ?/-links of the interval [{qx,l), {qx,%)) respectively. Note that = I 
when qx = I and Sy = I when qy = 1. Sg^, and 5*5,0 commute with Ly and all plaquette 
operators except the one located on the left of the origin (i.e. W^n^^i)) and Lx, with 
which they anti-commute. 

By using string operators, fermionic creation and annihilation operators are defined 
on each 2;-link q as 



4= 0, = ^^^^^^^ (8) 

It is not difficult to show that they satisfy the fermionic anti-commutation relations 

{Cj, Cg>} = 6g^g>, { C j , C^, } = {Cg , Cg>} = C (9) 

Because the c-fermions and the string operators have the same commutation and anti- 
commutation relations with the plaquette and the loop operators, the quadratic forms 
of fermionic operators ( e.g. CgC^, for any q and q') commute with Lx, Ly and Wg for all 
q. In other words, for a 2A^-spin system on a torus, the even fermionic states span the 
2^~^ dimensional subspaces of each {Wp, Lx, Ly} configuration. Therefore, only states 
having an even number of c-fermions are realized. 

By using these fermions, we write "H in Equation (5) as 

n=J2 J^^M - + C,+nJ (10) 

9 



+ Y,J];Yg{cl-Cg)icl^^^+Cg+, 

q 

+ Y^.r^{2clcg-i), 



Q 



where X and Y are defined (see Figure 2) as 



( ly-^ 

n W^fe,^.) ^/ ?x ^ Nx 

ty = l 

qy-1 

-Lx n W{q.,iy) if Qx = Nx 

iy = l 
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(b) (c) 

Figure 2. (Color online) (a) The transformation of a honeycomb lattice into a square 
lattice by contracting the ^-Hnks to a point, (b) and (c) show -^^(2,4), ^(3,4) respectively 
for a 4 X 4 square lattice whose opposite sites are identified. 



1 if Ny 

^'^''^^ ^ -Ly n n ^(....) ^/ qy = Ny. 



Similarly, the fermionic representation of the Pj^ terms of V (6) reads 

~ ~ '^q q+ny{.C^q+^y + Cg+n^ ) (c^_^-^_,_-__, + Cg+j^j^+fj^ ) , 
Pq = ~ ^^9(4 - Cg)(c|^_^_^^ - Cq+riy), 

^q ~ ~ '^q i'^g+rix{(^l+nx ~^ ^Q+'rix) {'^l+nx+riy ~^ ^Q+nx+riy), 

Pq = '^g iNq+nyyq+nxic^q-^i^y ~ ) (Cg+^^ — Cg+J^^ ) , 

Pq = Kq lXqYq{c\^^^ + Cg+fi J (cj_^^^ + Cq+Ux) ■ 

Let us define Pq and P^ as 

-Pg^ •= -Pg^ + Pq-ny = ~'^q ^2Xq (cjc^^jj^ + CgCg+ipiJ 
Pq '■= Pq + Pq-rix = ^'^g ^'^^Qi^q^l+ny + CgCg+n^) 

where we assume 1^^ = i^l = nj, 7; , i^n = i^l = f^t w • Now, we can write the V term as 

q q q tiy ' q q q^ii'X ' 

V=-Y.{P-q+P^q+Pl + Pt). (11) 

Notice that the total Hamiltonian T-itot = H + V of the Kitaev model is quadratic 
fermionic. The quadratic fermionic Hamiltonian H of any system with N fermions must 
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jk 



where ^ is Hermitian and A is antisymmetric, and can be rewritten in the following 
way [26] 



H 



(12) 



where 



c{ ...C]... cjv Ci ...C,... Cat 



ci ...Ci... cn 4 ^iv 



1 T 



Note that 



M 



-A* 



A 



(13) 



is a Hermitian matrix that can be diagonalized as M£, = T'^MT where T is a unitary 
operator of the form 



U V* 

V u* 



(14) 



whose columns correspond to eigenvectors of M. The matrix Md :- 



where 



E 
-E 

E is a diagonal matrix with positive entries. These are placed in an increasing order 
El < ... < i^AT as a convention [26]. By replacing M in Equation (12) with M = TMdT\ 
we get 



H 



1 

2 L 



4 



U V* 
V U* 



E 
-E 



Since T is a unitary matrix, it is possible to define new sets of fermions 



Pi 



U V* 

V u* 



(15) 



These two definitions are compatible, and therefore, the Hamiltonian H can be rewritten 
in a free fermionic form as follows: 



Ei 



(16) 
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For the total Hamiltonian T-Ltot of the Kitaev model, Aitot ^ the analog of the matrix 
(13) - is given in terms of X and Y; therefore it is diagonalized separately for each X 
and Y configuration (i.e. {Wp, L^, Ly} configuration). 

On the other hand, the eigenstates of the quadratic fermionic Hamiltonians can 
be explicitly written by using the Bloch- Messiah theorem [27,28]. According to this 
theorem, any unitary matrix of the form of T in Equation (14) can be decomposed into 
three matrices of very special form: 



(17) 



where D and C are unitary matrices and U and V are real matrices of the general 
block-diagonal form 



D 







U 


V 




' c 








D* 




V 


u 







c* 



u 











V 











where Z and / are the zero and identity matrices of the same size respectively, and 



Ui 
Ui 



Vi 



Vi 
-Vi 



where and Vi are positive real numbers. By using Equation (17) in Equation (15), we 
have 



4 



Here, D is used to define new operators a^" and a: 



D 







U 


V 




' c 








D* 




V 


u 







c* 



D 





H 







D* 




= c^D* 



Then there is a special Bogoliubov transformation 



u 


V 




V 


u 





This distinguishes the "paired" levels {up > 0; fp > 0) 



Oin 



V/pOp ~\~ VpOp^ 



Olp VpOp 



Upttp, 



(where {p,p) are defined by Ui and Vi (18)) from the "blocked" levels (w^ = 0; 



a' 



0!r. 
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which are either occupied {vi = 1; Ui = 0) or empty. Finally, a linear transformation of 
the and a by the unitary matrix C gives: 



Pi /3. 



' c 





-1 


f PI 


= ale 





c* 




[ p^ 


= a^C* 



For a general quadratic fermionic Hamiltonian, the ground state wavefunction is 
defined as a non-zero wavefunction |0) such that Pk\4') = for all k. It can be easily 
verified that the following wavefunction satisfies these criteria: 

i0)=n«Jn(^^'+^x4)i-)' (19) 

i p 

where |— ) represents the vacuum of c-fermions. In the honeycomb model, the vacuum 
|— ) belongs to the X,y-configurations for which Aitot is diagonalized. However, because 
only the even fermionic configurations are allowed for each X,y-configuration, when the 
number of elements in the first product (i.e. i-part) of (19) is odd, the ground state 
1$) of the model is the next excited state |$) = PI\4>) where PI is the minimum energy 
fermion. In that respect, the number of a] determines the fermionic parity of the system. 
For future reference, it is important to associate every eigenstate of the system with 
U V* 



a T 



V U* 



matrix such that the Bloch-Messiah representation Equation (19) 



represents that eigenstate. For an excited state this can be done by exchanging 

the roles of Pi{i) and Pl{i) which manifests itself as the exchange of the ith. and (A^ + z)th 
columns of T; we denote this matrix by T' . Thus the ground state \(f)') of T' is equal 
to A straightforward generalization of this method to the other excited states 

associates any eigenstate with a particular T matrix. 



3. Adiabatic Motion of Vortices and the non-Abelian Berry Phase 

In this Section, we discuss how to move vortices from one plaquette to another, then 
describe the evolution of the system under an exchange of vortices. 

By changing the sign of the relevant coupling coefficients J and k, we can emulate 
the fermionic spectrum relevant to any X,y-configuration starting from the trivial X, Y- 
configuration (i.e. Xq = 1 and Yq = 1 for all q) which we call the reference X, Y- 
configuration. 

This approach, first introduced in [18], allows us to consider the Hilbert spaces of 
different X, y-configurations connected by the coupling coefficients J and k. 

In this way, we can also simulate the creation, annihilation and motion of the 
vortices. For example, consider the {Wp, L^, Ly} configuration (see Figure 3(a)) 

= Ly = —1 and Wg = 1 for all q except Wq/ = Wq/^n^ = —1; 

where q' = (g^, q'y) such that q'y 1. This configuration gives Yg = 1 and Xq = 1 for all 
q except Xqi = — 1. The fermionic spectrum of this configuration can be achieved from 
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= (1,1) 



(a) 



= (1,1) 



(b) 



Figure 3. (Color online) (a) The creation of two vortices sharing an a;-link. (b) The 
creation of two vortices sharing a j/-link. 



the reference X, y-configuration by using the negative values of the set of coefficients 
[Jk]^> '■= {Jq'y i^qz-nyy '^q'y '^q'} (^66 (10) and (11)). lu othcr words, changing the sign of 
[Jk]^', := { Jg,, i^q'-riyy ^q'' ^q'} t)e considered as the creation of two vortices from the 
vacuum (i.e. the reference X, F-configuration with positive J and k values). 

For an analogous configuration shown in Figure 3(b), vortices on the plaquettes q' 
and q' — for q'y = Ny and q'^ ^ 1 can be created from the reference X, F-configuration 
by changing the sign of the set of coefficients [Jk]^, := {J^,, f^l'^^^, i^gi, i^l'}- 

Generally, changing the sign of [Jk]^ alters the vorticity of two plaquettes 
and when qy ^ 1, and may be seen as creation, annihilation or motion of 

vortices, depending on the initial vorticity of the plaquettes. A similar effect can be 
achieved on the plaquettes and for q satisfying qy = Ny and 7^ 1 

when we change the sign of [</k]^. We point out here for completeness that it is also 
possible to simulate the vortices by changing the [Jk]^ when g^^ = 1 or by changing 
[Jk]^ for q satisfying qy ^ Ny 01 q^ = 1 by carrying the sign to some of the c-fermions. 
However we will use the previous approach to exchange vortices as it is sufficient. 



1 
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o = (i,i) (iV.,1) 



Figure 4. (Color online) A configuration with 4 vortices having minimum distance 
d = 3 and sizes — 10, Ny = 6. Opposite sides of the lattices are identified. Red 
links highlight the adiabatically changed links. The path swapping the vortices B and 
C consists of the links on the arrows 1, 2 and 3. 
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From now on we will use the term vortex only for these simulated vortices. In Figure 
4 a template configuration is shown with four vortices which are created from the vacuum 
by changing the sign of [JnY x-links between A and B, and between C and D 

(colored by red in Figure 4). In this paper, we will work with several configurations 
similar to Figure 4 with different sizes as: A^^^ = 3d + 1, Ny = 2d when d is odd 
and Nx = 3d, Ny = 2d when d is even for d = 1,...,9, where d is the minimum 
distance between the vortices. Note that all configurations are even-by-even; the other 
configurations (odd-by-odd, even-by-odd etc.) will be studied in future. For all these 
configurations, the coupling coefficients for the vacuum are identical for all plaquettes: 
= J| = J| = J = 1 for all q. The strength of the effective magnetic field k for the 
vacuum is also the same for all the plaquettes q, = k,"^ = k,^ = = k, and is taken 
from two sets which differ in magnitude: (i) k = / x 0.01 where / = 1,...,5, and (ii) 
K = I X 0.05 where / = 2, 8. 

To swap the position of the vortices B and C - see Figure 4 - we need to adiabatically 
move the vortices along the paths indicated by the arrows 1, 2, and 3. This requires 
us to slowly change the sign [Jk] of the links which are intersected by the paths. For 
the large values of k, i.e. the values from the set (ii) above, all these configurations 
have a (nearly) two-fold degenerate ground state which is separated by a gap from the 
rest of the Hamiltonian's spectrum. This can be seen from Figure 5(a) which shows the 
average of the ratio of the splitting of the ground state degeneracy and of the gap to the 
first excited state G'2,i(A)/G'3,2(A) along the path. Here G„,^(A) = |E"^(A) - ^"(A)| is 
the energy difference between the m*'* and n*'' eigenstates and A parametrizes the path, 
i.e. it represents the values of [Jk] for the links on the path. The degeneracy decreases 
with the minimum distance d. The ratio G'2,i(A)/G'3.2(A) for small values of n (i.e. from 
the set (i)) exhibits more involved behavior as the gap between the ground states and 
higher excited levels is much smaller in this case. 

It is interesting to point out that fermionic parity of the system may change while 
the particles are exchanged. The average values of the parity of the system are shown 
in Figure 5(b) where 2 and 1 are assigned to the even and the odd parity sectors 
respectively. The parity of the system sets the small gap (72,1 between the nearly 
degenerate states as ^2,1 = — Ei for odd systems and 6*2,1 = E2 + Ei for even 
systems, where Ei denotes the spectrum of the system increasing with the index i > 1. 
However, the parity does not affect Gm.2 for any value of m. This is one of the causes 
of the oscillations seen in Figure 5. It is not the only cause though. Calculations show 
that the average energy of the nearly-zero modes {Ei and E2) oscillates in the same way 
as well. 

For a system which is initially in the ground state space, the exchange of B and 
C transforms the system within the 2-dimensional ground state space provided that 
the process is adiabatic (Appendix A contains a detailed discussion on adiabaticity and 
the path which is followed). When the change of the parameters of the Hamiltonian 
follow a smooth closed curve in the parameter space, then the evolution of the system 




Figure 5. (Color online) (a) Average value of the ratio G2,i(A)/G3,2(A) along the 
trajectory vs. minimum distance d between vortices, (b) The average values of the 
parity of the system while the particles are exchanged vs. d. Average oddness/evenness 
of the systems along the trajectory is calculated after assigning 2 and 1 to the even and 
odd parity systems respectively. Note that the numerical details of the calculations 
are given at the beginning of Section 3. 
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is governed by the following Berry phase (matrix) B{C) [30]: 

B{C) := Pexp |z ^A{X)dX^ (20) 
:= lim exp {z^(Aa/)AA} ■ ■ -exp {z^(Ai)AA} 

M-^oo 

where V denotes the path-ordered integral, Aabi^k) = i Ia-a^' — 

{1,2}, and Ai and Aj\/+i are coinciding points denoting the beginning and end point of 
the closed trajectory whose curve length length(AA/+i, Ai) is divided into M equal pieces 
AA =length(AM+i, Ai)/M. Note that because the path traveled by vortices encloses zero 
area, all local and geometrical effects are eliminated and only the topological interaction 
is realized. 

The Berry matrix B{C) is written in the same basis as that of ^(Ai) (see 
Appendix A for more details about Ai). However, it is independent of the choice of 
the bases in which other ^(A)s are written, as long as the basis states are smooth 
functions of the parameter A. For practical purposes, we used the nearly degenerate 
energy eigenstates ([^^(A)) and |$^(A))) of the system. We point out for clarity that 
|$^(A)) and |$^(A) are always distinguishable thanks to the small gap which separates 
them and which is never less than 10~^ for all values of k and d we used. 



4. Numerical Methods for Calculating the non-Abelian Berry Phase 



In this Section, we present the arguments that we use for the numerical calculation of 
the Berry matrix. The results and discussions are presented in the next Section. 
In order to calculate the Berry matrix i3(C) (20), we first need to evaluate 



Aab{\k) = I mx)\j^^\x)) 



a,6 = {l,2}. 



(21) 



There are many different ways to approximate the derivative of a function [31], but we 
are going to use the central-difference formula which says 

/(x + Ax) -/(x- Ax) 



/'(x)Ax 



0((Ax) 



(22) 



as long as the third derivative of / is continuous. It approximates /'(x) on the order 
0((Ax)^); however it is also possible to get higher order approximations [31]. In our 
case, we performed the calculations of /'(x) on the order 0((Ax)^). 
By applying the formula (22) to Aab{\), we get 



■^abi^k) 



('l'"(A,)|'l'^(Afc+i)) - ('l'"(Afc)|$^(Afc_i)) 
2AA 



This reduces the Berry matrix calculation to finding the overlaps between the ground 
states of adjacent points on the trajectory. 

Let M be the number of data points used to calculate the Berry matrix (20). For 
all data points k = 1, M and some other point /c = 0, let be the ground states 
of /3(/c)-fermions defined as 



T{k) 



(23) 
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where T{k) := 
we can write 



U{k) V*{k) 
V{k) U*{k) 



for A; = and replace 
(24) as 

^ Pl{k) P^{k) 



as in equation (14). Since the T's are unitary matrices, 

/3t(0) /3^(0)]rt(0) (24) 
of Equation (23) with the left hand side of Equation 



where 



and 



T(A;,0) :=Tt(0)T(A;) 



/3t(0) /3«(0) T(A;,0) 



f/(A;,0) \/*(A;,0) 
\/(A;,0) f/*(A;,0) 



(25) 



(26) 



U{k, 0) := f/^(0) U{k) + Vt(0) r(A;) 
\/(A;, 0) := V^iS)) U{k) + U^{0) V{k). 

By applying the Bloch- Messiah theorem to Equation (25), we can only get the 
absolute value of the overlap between the ground states |(i;Ao|0A;)| = \/| det U {k, 0)| (see 
Appendix B). However, it is possible to get the complete overlap by using the Thouless 
theorem [28,32]. The Thouless theorem states that when (^fcl^o) 7^ 0, the ground state 

can be written as 

\<Pk) = |^(fc,o)) := V|dett/(A:,O)|e^(^^'°)|0o) (27) 

where 

Z(^,0) = ^5^Z„„,(A:,0)/3t(0)/3i,(0) and Zik,0) = {Vik,0)U~\k,0)y . 

n,n' 

Having the ground state represented as \ip{k,o)) for k = 1, M where we have 
fixed the overall phase to 1 (Equation (27)), we recall that the excited states can be 
represented by using the column exchange technique on T{k, 0) discussed at the end of 
Section 2. Then, the overlap between the ground states \ip(ifl)) and \ip{kfi)) for /c 7^ / 7^ 0, 
reads 

{i^mli^im) = VI detf/(/,0)| V| detf/(A;,0)| (0o|e^'(''°)e^('='°)|0o). 
In a recent paper [33], the overlap (0o|e^^(^'°)e^('='°)|0o) has been calculated as 

Pf(Z(/, 
Z{k,0) 



where Pf denotes the Pfaffian, Z{1, 0; /c, 0) 



and N is the 



-/ 

/ -Z*{1,0) 

number of fermions (and therefore also the size of Z). Since the numerical algorithms 
for calculating the Pfaffian of large matrices are generally slow, it is more convenient to 
proceed using the following expression (see Appendix C) 

= (-l)^(^+i)/2 v/exp{2^o(/,fc)} \ detUil,k)\ (28) 
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where 

U{1, k) := U\k)U{l) + V\k)V{l) 

eo{l,k) = aTg{detU{k,0)detU^{l,0)detU{l,k)} . 

Although the sign in Equation (28) is ambiguous due to the square root, for a series of 
smoothly varying matrices the correct sign can be traced from their previous values. 

We point out here that, for the purpose of calculating the Berry matrix (20), the 
degenerate ground states at each point of the trajectory in parameter space is defined 
in terms of the reference state |0o)- Thus the reference state must be chosen in such a 
way that it is not orthogonal to any of the ground states along the whole trajectory, as 
in Equation (27). 

5. Numerical Results and Discussions 

Before presenting the numerical results, let us briefly discuss the Ising anyons consisting 
of three different particles: 1 (i.e. vacuum), e and a. Bringing two particles together is 
called fusion, and Ising anyons satisfy the following fusion rules: 

a X a = 1 + e, a x e = a, exe = l, (29) 
1x1 = 1, 1 X a = a, 1 X e = e 

where x denotes the fusion operator. For example, the first rule says that two a- 
particle may either annihilate or fuse into an e-particle. Although the fusion of two 
Abelian anyons gives one outcome, non-Abelian anyons have multiple fusion possibilities 
or channels which are one way of accounting for the degeneracy of the ground state. 
Thus (T-particles are non-Abelian anyons while the other two particles are Abelian. 

To detail the nontrivial implications of these rules, consider a system with four 
(T-particles (a, b, c and d) whose total topological charge corresponds to the vacuum. In 
this system, the fusion of any two cr-particles (say a and b) determines the fusion channel 
of the other two cr-particles (c and d) and because there are two different fusion results 
that a and b can fuse into, there is a two dimensional fusion space associated with the 
system. The basis of this space can be chosen as the resulting fusion states of a and b 
as On the other hand, another basis can be chosen based on the fusion 

results of b and c {Ifi^"^), The matrix that relates the two fusion bases is called 

the F-matrix. It is analogous to the 6j symbols encountered in the couplings of three 
spin-1/2 particles. For arbitrary numbers of anyons, different sequences of fusion-basis 
transformations, starting from a particular fusion-basis ending in another one, must 
be equal. This imposes a consistency condition on F-matrices known as a pentagon 
equation [11]. 

On the other hand, particle exchange operators r are represented by /^-matrices 
on this fusion space. Nontrivial relations arise when we consider particle exchange 
operators together with F-matrices. These relations can be expressed by the hexagon 
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equations [11]. Solving the pentagon and hexagon equations specify the consistent anyon 
models. For the Ising model, there are consistent solutions as follows 

i?f = -1, i?^' = R'f = la, i?r = ^e^""/^ R^/ = Qe-'^'^l'^ 

with the possible combinations of topological spin Q (associated with the counterclock- 
wise rotation of a particle by angle In around itself) and a 

e = a = (-i)(-+i)/2 

where v is the spectral Chern number taking odd integer values [11]. For translationally 
invariant configurations of the Kitaev Honeycomb model (e.g. the reference X, Y- 
configuration i.e. Xq = 1 and Yg = 1 for all g), the spectral Chern number is equal 
to ±1 (depending on the direction of the magnetic field) [11]. 

The solution of the pentagon equation also determines the transformation relation 
between the fusion bases If^"'')} ^^"^ { I ^1^)5 l^e"^)} of cr-particles (a, 6, c and d) as 

up to an arbitrary relative phase e*"^. 

In that regard, the representation TZ of the braiding operator which exchanges 
cr-particles h and c is diagonal for the fusion basis 

where {he} stands for the basis dfi^'^), l^^e*^)} and by using the transformation relations 
above, 7^{a6} reads 

^ 1 / i?r + R7 e-n^r - R7) \ .3. ^ 

{ah} 2 \^ e'^{Rl'' - R^'') Rl'' + RI'' j' ^ ^ 

On the other hand, the Berry matrix i3(C) given by Equation (20) is written in a 
different basis. This is the same basis as ^(Ai), where Ai is the starting point of the 
trajectory C in the parameter space. These basis states consist of two nearly degenerate 
energy eigenstates {|<l>^(Ai)), |$^(Ai))} (see Sec. 3 also). We mention for completeness 
that the coupling coefficient of the starting configuration is slightly different than that 
of the original configuration shown in Figure 4. We discuss the details of the trajectory 
used in the calculations in Appendix A. 

Because the Berry matrix is not invariant under the basis transformation and the 
relation between the fusion channels and the basis {|$^(Ai)), |$^(Ai))} is unknown, the 
comparison of the eigenvalues of the Berry matrix with {Rl" ■, R^} is more meaningful. 
However, before we do that we would like to point out an interesting similarity between 
the calculated Berry matrices and an 7^-matrix: by making an analogy between four 
cr's (a, 6, c and 0?) and four vortices (A,B,C and D) in Figure 4, for the values of n from 
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Figure 6. (color online) |i3(C) — 7?.{ab}|p : Frobenius distance between B{C) and 
d. Note that the numerical details of the calculations are given at the 



n 



{AB} 



vs. 



beginning of the Section 3. 



the set (ii), the numerical Berry matrix B{C) converges to TZ{ab} with spectral Chern 
number z/ = 1, namely, 

g— 7i7r/8 givr/f 



^{ab} = 




as the minimum distances d increases. We point out that the arbitrary relative phase 
e**^ = —i is chosen for TZ^^}', and a similar arbitrary relative phase is fixed between the 
basis states |$^(Ai)) and |$^(Ai)). This arbitrary phase refiects a gauge freedom and is 
chosen to obtain the best agreement between B{C) to 7^{ab}- 

The numerical results are summarized in Figure 6 where the Frobenius norm was 
used to measure the distance between the two matrices. The Frobenius norm of an 



1/2 



n X n matrix A is defined as \A\f = [^ij l^fcl^j = v^tr(y4y4t). Here, we calculate 
the norm |7^{ab} — -^(0)1^ to measure the distance between B{C) and 7^{ab}- The high 
level precision of the results is refiected in the unitarity of the calculated Berry matrix: 
|i3(C)i3^(C) — /|^ < 10~^ for all d and k. It is not difficult to show that the maximum 
Frobenius distance between two unitary 2x2 matrices is equal to 2, so to evaluate the 
error of the calculations the Frobenius distance has to be divided by 2. 

As an example, for the vortex configuration illustrated in Figure 4 with d = 9 and 
K = 0.25, the numerical value of the Berry matrix B{C) is 



BiC) 



0.653270 ^ 
-0.653280 



0.270630i 
- 0.270598i 



0.653280 + 0.270598i 
0.653296 + 0.270568i 



whose Frobenius distance to 7^{ab} given in Equation (5) is equal to 4.8 x 10 ^ and the 
normalized error associated with B{C) compared to the exact result from the effective 
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field theory 7^{^B} is 0.0024%. 

This similarity between the basis {|$^(Ai)), |$^(Ai))} and jl]^^)} can be 

understood from the perspective of Majorana fermions. First of all, note that the 
eigenstates 1$^) and |$^) are also the eigenstates of the operators i'jia'Jib and i'j2al2b 
where 

7ia = /3l + /3i, 7i6 = ^(/3l-/3i), 

are Majorana fermions, which are their own anti-particles. Kitaev showed that for odd 
values of the Chern number vortices carry Majorana fermions [11]. The localization 
of the Majorana fermions around vortices are also demonstrated numerically in [34] by 
expressing them in terms of c-fermions. In that respect, the equality of the bases can be 
understood as the Majorana fermion pairs (71a, 71;,) and (72a, 726) being localized on the 
vortex pairs (A, B) and (C, D). Note that because the localization of Majorana fermions 
is only related to the eigenstates of the Hamiltonian, it does not depend on the history 
of the vortex motion. That is, the configuration properties (such as relative distance of 
vortices to each other, boundary conditions, the values of and Ly etc.) determine 
the relation between the energy basis and the fusion basis. For example, for the same 
size configurations with = 1 and Ly = —1, the Berry phase B{C) is diagonal so that 
{|$i(Ai)), |$2(Ai))} ^ {|fif ), l^f )}• (Berry phase calculations for the other L^, Ly and 
N^, Ny configurations will be presented in a different context in future.) 

On the other hand, the eigenvalue comparison is presented in Figure 7 where BoiC) 
is the diagonalized form of B{C). To emphasize the proximity of the eigenvalues of B{C) 
(l.h.s.) to those of 7^{ab} (r.h.s.), the explicit values for d = 9 and k = 0.30 are 

0.3826813 + 0.9238804Z ~ e^*^/^ 
0.9238802 - 0.3826817i ~ e~*^/^ 

giving the Frobenius distance 3 x 10~^ and the error 0.00015%. 

The exponential convergence of the Berry phase to the expected values as 
the minimum distance d between the vortices increases, which the accuracy of our 
calculations allows us to observe, reveals subtle connections. The oscillations and the 
agreement of the calculated Berry phase to the expected values are in correlation with 
the exact ground state degeneracy and the gapfulness of the system (Figure 5) with 
few exceptions. However, one should keep in mind that although the ground state 
degeneracy and the gapfulness of the system play a role in the physical evolution of the 
system, the Berry phase only depends on the eigenstates. Therefore, all the Berry phase 
values should be understood in terms of the properties of the eigenstates. 

At this point, it is essential to note that the structure of the eigenstates determines 
the spread of the localization of the Majorana fermions around the vortices [34]. Because 
the braiding properties of the non-Abelian anyons are determined when the particles are 
away from each other to avoid tunneling, the radius of the Majorana fermions around the 
vortices affects the calculated Berry matrices. Thus, the relation between the calculated 
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Figure 7. (color online) |i3D(C) — 7?.{bc}|j;^ • Frobcnius distance between diagonalized 
Bd{C) and TZ^^c} vs. d. Note that the numerical details of the calculations are given 
at the beginning of the Section 3. 



Berry phase to the expected TZ values is also an indication of how tightly the Majorana 
fermions are localized around the vortices. 

In that respect, the correlations between Figure 5 and Figure 7 suggests that the 
bigger gaps and the nearer degenerate states make the localization of the Majorana 
fermions tighter around the vortices. However, it seems that this is not always the case 
as is seen for the systems with d = 6 and the small values of k (i.e. from the set (i)). 

6. Outlook and Summary 

To summarize, we have analyzed some four-vortex configurations of the Kitaev 
honeycomb model (see Figure 4) under various magnetic fields. All lattice configurations 
are even-by-even in terms of the number of plaquettes and they are defined on a torus 
whose size is determined by the minimum distance d between the vortices. As was shown 
in Section 3, the ground states of such configurations are two-fold nearly degenerate. 
Under the adiabatic interchange of the vortices, the evolution of the system is restricted 
to the ground states and it is governed by the Berry matrix whose numerical evaluation 
was presented in detail in Section 4. The numerical results show that as d increases the 
small gap between nearly degenerate state decreases exponentially with some oscillations 
(see Figure 5) also the calculated Berry matrices get exponentially closer to the expected 
ones with similar oscillations (see Figure 7). 

The main result of the paper is the explicit demonstration of the non-Abelian 
statistics by numerically calculating the Berry matrix which governs the evolution 
of the system under the adiabatic interchange of vortices of the Kitaev honeycomb 
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model. We showed that the resulting Berry matrices exponentially converge to the 
expected braiding properties of Ising anyons derived from the effective field theory. 
The presented method for the Berry phase calculations represents a direct approach 
to the non-Abelian statistics, and it can be used to study braiding properties of 
anyons for any configurations. The high accuracy of the method also allows us to 
test the stability of the Berry phase, so that we may make conclusions about both its 
implementation under (more) realistic conditions and the pitfalls we should expect in 
such situations. Moreover, the accuracy of the calculations presented here can provide 
important predictive power in the analysis of other, less understood, topological phases. 
They also reveal the dependence on details of the implementation of the braiding 
operations and thus allow testing of non-adiabatic and other effects which will likely 
constitute a dominant source of errors in topological quantum information processing. 
The accuracy of the calculations shows the exact dependence of the simulated Berry 
phase on the details of the model, like the splitting of the ground state levels intrinsic 
to any finite system. Possible applications thus extend to modeling and implementation 
of quantum information protocols whose reliability can be tested under various effects; 
for example, disorder. Naturally, the first step in this potentially fruitful story is the 
demonstration of highly accurate and sufficiently large scale direct simulation of the 
Berry phase as we have presented in this manuscript. Moreover, the method is general 
enough to be useful in the study of non-Abelian statistics in other models including 
the Yao-Kivelson model [21], square-octagon model [22] or any system with a quadratic 
fermionic Hamiltonian. 
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Appendix A. 

When we change the parameters of the Hamiltonian so as to follow a smooth closed 
curve in the parameter space, the evolution of the system is governed by the Berry 
phase if the process is adiabatic [30]. In this Appendix, we will introduce the trajectory 
which we use and then discuss the adiabaticity of the exchange process in detail. 

To examine the smoothness of the trajectory shown in Figure 4, let us take as an 
example the [Jk] of the first two links of the path, namely [Jk]^ and [Jk]^^- where 
q = {2,Ny). We need to change [Jk]^ from a to —a to move vortex B to the right 
plaquette where a = 1 for J, and a = I x 0.01 or a = k x 0.05 for n where / = 1, 5, 
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(a) 



(b) 



Figure Al. (Color online) (a) Changing the J and k values of a link after another 
gives us a trajectory with square edges, (b) Changing the J and k values of the next 
link before stopping to change that of the present one gives us a smooth trajectory. 



k = 2, ...,8. Changing the sign of [Jn]g first, and then [J^j^^j^^, results in a trajectory 
with discontinuous edges as shown in Figure Al(a). To move the vortex smoothly, we 
need to start changing [Jk]^_^_^ before [Jn]g reaches —a, as shown in Figure Al(b). For 
this reason a vortex is never completely localized on a particular plaquette. We will refer 
to the part where only one [Jk] changes as the linear part and the part where two [Jk] 
changes as the circular part. In our numerical calculations, every link on the trajectory 
is sampled with 4000 data points (3991 linear plus 9 circular) which are separated by an 
equal distance As in the parameter space. Let / be the length of the linear part from 
to the beginning of the circular part (see Figure Al(b)) and r the radius of the circular 
part; then r + I = \a\ and As = Al = rAd where Al = //3991 and A9 = f /9, so that 
/ = 0.997|a| and r = 0.003|a|. 

The starting point of the trajectory (i.e. Ai) is chosen to be the beginning of the 



linear part of the [Jk] of the first link (i.e. [Jk 



\y 



of the path; because of that the 



coupling coefficients of the starting configuration are slightly different than that of the 
original configuration shown in Figure 4. 

On the other hand, the first requirement for a process to be adiabatic is the gap 
between the eigenspace and the rest of the Hamiltonian's spectrum and this has been 
analyzed in Figure 5. The standard condition for maintaining the adiabaticity for the 
j^th eigenstate is as follows: 



($"1$™) 



En - 



<1, 



where dot denotes the time derivative. However, this condition is only valid if 
($"'(A)|$™(A)) and £'™(A) - i?"(A) are constant for all m through the trajectory [29]. 
When ($"(A(t))|<i>'"(A(t))) and E™(A) - E"(A) are not constant the validity condition 
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Figure A2. (Color online) Adiabaticity condition minG3 2(A) where 

max|^$(i'2)(;j^)|^ := max(|^$i(A)|_L , ^ The numerical details 

of the calculations are given at the beginning of the Section 3. 



of the adiabatic approximation can be written as 
max(E^^J(A$n(A)|$-(A))| 



min|E"(A) -E™(A)| 



< 1. 



where t; = ^ is the speed of the process and we used the following equality 
($"(A(t))|<i>™(A(t))) = -($"(A(t))|$-(A(t))) 

which can be verified by differentiating the orthogonality condition ($"'(A)|$"*(A)) = 0. 
First we define 

^|(A$".(A)|$«.(A))| (A.l) 
i=i 

for a /c-fold degenerate eigenspace with basis |$"'(A)) for i = 1, ...,k. Then by using the 
following inequality 



A$"(A)|$-(A)) 



< 



we can rewrite the general validity condition for adiabaticity as 

< 1. 



max 


A$-(A) 




min 


^"(A) - 


(A)| 



Figure A2 shows the validity of the adiabatic approximation for various k values 
against d. The smaller the values of k, the smaller the gap to the first excited state, 
and thus the smaller the speed v is required to maintain adiabaticity. 
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Appendix B. 



Here, we would like to show by using the Bloch-Messiah theorem [28] that the absolute 
value of the overlap K^ol^fc)! between the ground states of /3(0) and /3(/c)-fermions can 
be calculated as |(0o|0fc)| = a/| det U {k, 0)|. Let us start from Equation (25): 



where 



T{k,0) ■=T\0)T{k) 



(31(0) (3^(0) T{k,0) 



U{k,0) V*{k,0) 
V{k,0) U*{k,0) 



(B.l) 



and 



U{k, 0) := f/^(0) U{k) + 0(0) V{k) 
V{k, 0) := V^{0) U{k) + U^{0) V{k). 
The Bloch-Messiah decomposition of T{k, 0) reads 



T{k,0) 



D{k,0) 
D*{k,0) 



U{k,0) V{k,0) 
vlk,0) U{k,0) 



C{k,0) 
C*{k,0) 



Let us define /3(0)-fermions as 

t{0) = Y,Dk'k{k,0)/3l{0) 

k> 

Now, we can write the vacuum |0(fc,o)) of /3(/c)-fermions in terms of the vacuum |0o) of 
/3(0)-fermions as 

= n^Jw n k(fc,o) +^,(fc,o)^i:(o)^i(o)) i^o). (b.2) 



Note that and |0(A:,o)) are the same up to an overall phase. 
Now the overlap (0o|0{fc,o)) reads 

(0o|0(M)) = (0o| WJ\ n (np(A;,0) + ^;,(A;,0)^];(0)4(0)^ 

i p 

This is only non-zero if there is no product with index i and it is equal to 

(0o|0(M)) = (0oin^p(^'O) l<^o) = n^p(^'O) = \/detf7(/c,0). 

p p 

Recall that U{k,0) is diagonal with positive entries. Because of the arbitrariness of 
the phase of Equation (B.2), the overlap {(po\4>{k,o)) can be taken to be positive; hence 

|(0o|0fc)| = (0o|0(fc,O))- _ 

Moreover, since U{k,0) = D{k,0)U{k,0)C{k,0), and D{k,0) and D{k,0) are 



unitary, ^JdetU{k,0) = ^/\ detU{k,0)\. Therefore, 
(0o|0fc) = A/|detf/(A;,0)|. 



(C.l) 
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Here, we would like to show that the overlap 

{i^m |^(fc,o)) = V|dett/(/,0)| V|det[/(A:,0)| pf(2(/, Q; A;, 0)) 

where 

^ Z{k,0) -I 
I -Z*{1,0) 

and is the number of fermions and therefore also the size of Z matrices, can be written 

as 

ii^mlAm) = Vexp{z^o(/,A;)} \ detU{l,k)\ 

where 

U{1, k) := U\k)U{l) + V\k)V{l) 

Ooil.k) = 8.Tg{detU{k,0)detU\l,0)detU{l,k)} . 

First of all note that for any skew-symmetric matrix K, Fi(K) = ^ydet{K). 
~ A B~ 
C D 

coefficients and CD = DC, then detfsT = det {AD - BC) [35]. Therefore, 



Moreover, if K 



where A,B,C,D are n x n matrices with complex 



Pf(Z(/, 0; k, 0)) = ^/det Z{1, 0; k, 0) = ^det (/ - Z{k, 0)Z*{1, 0)). (C.2) 

Recall that Z{k, 0) = {V{k, 0)U~^{k, 0))*. By using the conditions that U{k, 0) and 
V{k, 0) need to satisfy in order for the matrix T{k, 0) to be unitary, it can be shown 
that Z{k, 0) is also skew- symmetric and can be written as 

Z{k,0) = - {u\k,o)yW\k,o). 

Therefore, Equation (C.2) reads 

/ - Z{k, 0)Z*{1, 0) = / + {U\k, 0))-^V\k, 0)V{1, 0)U'\l, 0) 

= {u\k, o)y^ {u\k, o)u{i, 0) + v\k, o)v{i, o)) u^\i, o). 

By using Equation (26), it is easy to show 

U\k, 0)U{1, 0) + V\k, 0)V{1, 0) = U {I, k) , 

where U{1, k) := W{k)U{l) + V'^{k)V{l) is a generalized version of Equation (26). 
Therefore we have 

det (J - Z{k, 0)Z*{1, 0)) = det {{U^ik, 0)y^U{l, k)U~\l, 0)) . 

Finally, we can express {(f)Q\e'^''^'''^^e'^^^'^^\(f)Q) as 



zt(/,o)^z(fc,o)u \ — ( n^(^'+i)/2 / det U{1, k) 

' ^/detf/t(A:,0)detf/(/,0)^ 
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and express the overlap between ^^(kfi)) and o)) as 



(^a,o)l^(fc,o)) 



v/|dett/(/,0)| V|dett/(A:,0)| 



det W{k,0) detU{l,0) 



det U{1, k) 



^_^)7V(7v+i)/2 ^exp{zeo(/,A;)} |det[/(j, z) 



where 



^o(/,fc) = arg{det?7(A;,0) det 0) det fc)} . 
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